# 加载必要的包
library(tidyr)
library(randomForest)
library(ggplot2)
library(caret)
library(pROC)
library(gridExtra)
library(scales)

# ========== 1. 数据准备（假设Ttrain和Ttest已存在） ==========
# 请确保你的训练集数据框名为Ttrain，测试集数据框名为Ttest
# 如果名称不同，请修改下面的变量名

# 将目标变量转换为因子类型（训练集）
Ttrain$`营养不良风险` <- as.factor(Ttrain$`营养不良风险`)

# 将测试集的目标变量也转换为因子（用于后续评估）
Ttest$`营养不良风险` <- as.factor(Ttest$`营养不良风险`)

# ========== 2. 构建随机森林模型 ==========
# 使用训练集训练模型，ntree可根据需要调整
rfcla <- randomForest(`营养不良风险` ~ ., data = Ttrain, ntree = 200, proximity = TRUE)

# 查看模型摘要
print(rfcla)

# ========== 3. 可视化训练误差随树数量变化 ==========
# 提取误差数据并可视化
trainerror <- as.data.frame(rfcla$err.rate)  # 直接提取误差矩阵
colnames(trainerror) <- paste("error", colnames(trainerror), sep = "")
trainerror$ntree <- 1:nrow(trainerror)
trainerror <- gather(trainerror, key = "Type", value = "Error", -ntree)  # 修正：排除ntree列

# 绘制误差变化图
ggplot(trainerror, aes(x = ntree, y = Error)) +
  geom_line(aes(linetype = Type, colour = Type)) +
  theme(legend.position = "bottom") +
  ggtitle("随机森林分类模型 - 营养不良风险预测") +
  theme(plot.title = element_text(hjust = 0.5))

# ========== 4. 变量重要性分析 ==========
# 查看变量重要性图
varImpPlot(rfcla, main = "营养不良风险预测 - 变量重要性图")

# 输出变量重要性排序
importance_scores <- importance(rfcla)
print("变量重要性排序（基尼系数下降）：")
print(importance_scores[order(importance_scores[,1], decreasing = TRUE), , drop = FALSE])

# 可视化变量重要性（Top 15）
importance_df <- as.data.frame(importance(rfcla))
importance_df$Variable <- rownames(importance_df)
importance_df <- importance_df[order(importance_df$MeanDecreaseGini, decreasing = TRUE), ]

ggplot(head(importance_df, 15), aes(x = reorder(Variable, MeanDecreaseGini), y = MeanDecreaseGini)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  coord_flip() +
  labs(title = "Top 15 重要变量 - 营养不良风险预测",
       x = "变量", 
       y = "平均基尼系数减少量") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

# ========== 5. 测试集预测 ==========
# 预测类别
Ttestpre <- predict(rfcla, Ttest)

# 预测概率（用于ROC曲线）
Ttestprob <- predict(rfcla, Ttest, type = "prob")

# 查看预测结果频数分布
table(Ttestpre)

# ========== 6. 模型评价指标 ==========
# 确保预测标签和真实标签的因子水平一致
Ttestpre <- factor(Ttestpre, levels = levels(Ttest$`营养不良风险`))

# 混淆矩阵和各项指标
conf_matrix <- confusionMatrix(Ttestpre, Ttest$`营养不良风险`)

cat("\n========== 混淆矩阵 ==========\n")
print(conf_matrix$table)

cat("\n========== 模型评价指标 ==========\n")
cat(paste("准确率 (Accuracy):", round(conf_matrix$overall['Accuracy'], 4), "\n"))
cat(paste("95% 置信区间:", round(conf_matrix$overall['AccuracyLower'], 4), 
          "-", round(conf_matrix$overall['AccuracyUpper'], 4), "\n"))
cat(paste("Kappa值:", round(conf_matrix$overall['Kappa'], 4), "\n"))
cat(paste("灵敏度 (Sensitivity/Recall):", round(conf_matrix$byClass['Sensitivity'], 4), "\n"))
cat(paste("特异度 (Specificity):", round(conf_matrix$byClass['Specificity'], 4), "\n"))
cat(paste("阳性预测值 (Precision):", round(conf_matrix$byClass['Pos Pred Value'], 4), "\n"))
cat(paste("阴性预测值 (Neg Pred Value):", round(conf_matrix$byClass['Neg Pred Value'], 4), "\n"))
cat(paste("F1分数:", round(conf_matrix$byClass['F1'], 4), "\n"))

# ========== 7. ROC曲线与AUC ==========
# 创建ROC对象（使用第二列概率，即预测为“1”的概率）
roc_obj <- roc(Ttest$`营养不良风险`, Ttestprob[, 2])

# 计算AUC及其置信区间
auc_value <- auc(roc_obj)
auc_ci <- ci.auc(roc_obj)

cat(paste("\nAUC值:", round(auc_value, 4), "\n"))
cat(paste("AUC 95% CI:", round(auc_ci[1], 4), "-", round(auc_ci[3], 4), "\n"))

# 获取最佳阈值（约登指数）
optimal_idx <- which.max(roc_obj$sensitivities + roc_obj$specificities - 1)
optimal_threshold <- roc_obj$thresholds[optimal_idx]
optimal_sensitivity <- roc_obj$sensitivities[optimal_idx]
optimal_specificity <- roc_obj$specificities[optimal_idx]

cat(paste("最佳阈值:", round(optimal_threshold, 4), "\n"))
cat(paste("最佳阈值下的灵敏度:", round(optimal_sensitivity, 4), "\n"))
cat(paste("最佳阈值下的特异度:", round(optimal_specificity, 4), "\n"))

# 准备ROC绘图数据
roc_data <- data.frame(
  fpr = 1 - roc_obj$specificities,
  tpr = roc_obj$sensitivities,
  threshold = roc_obj$thresholds
)

optimal_point <- data.frame(
  fpr = 1 - optimal_specificity,
  tpr = optimal_sensitivity,
  threshold = optimal_threshold
)

# 绘制ROC曲线
p_roc <- ggplot(roc_data, aes(x = fpr, y = tpr)) +
  geom_line(color = "#2C7BB6", size = 1.5) +
  geom_abline(linetype = "dashed", color = "gray50", size = 0.8) +
  geom_point(data = optimal_point, aes(x = fpr, y = tpr),
             color = "#D7191C", size = 5, shape = 19) +
  scale_x_continuous(breaks = seq(0, 1, 0.2), labels = percent,
                     minor_breaks = seq(0, 1, 0.1)) +
  scale_y_continuous(breaks = seq(0, 1, 0.2), labels = percent,
                     minor_breaks = seq(0, 1, 0.1)) +
  labs(
    title = "随机森林模型ROC曲线",
    subtitle = paste0("营养不良风险预测 | AUC = ", round(auc_value, 4),
                      " [95% CI: ", round(auc_ci[1], 3), "-", round(auc_ci[3], 3), "]"),
    x = "假阳性率 (1 - 特异度)",
    y = "真阳性率 (灵敏度)",
    caption = paste("最佳阈值:", round(optimal_threshold, 3),
                    "| 灵敏度:", round(optimal_sensitivity, 3),
                    "| 特异度:", round(optimal_specificity, 3))
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 11, color = "gray30"),
    plot.caption = element_text(size = 9, color = "gray40", hjust = 0.5),
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.grid.major = element_line(color = "gray90"),
    panel.grid.minor = element_line(color = "gray95"),
    panel.border = element_rect(fill = NA, color = "gray70", size = 0.5),
    coord_equal()
  ) +
  xlim(0, 1) + ylim(0, 1) +
  coord_equal()

# 显示ROC图
print(p_roc)

# ========== 8. 预测概率分布图 ==========
prob_data <- data.frame(
  True_Class = Ttest$`营养不良风险`,
  Probability_1 = Ttestprob[, 2]
)

p_prob <- ggplot(prob_data, aes(x = Probability_1, fill = True_Class)) +
  geom_histogram(alpha = 0.7, bins = 30, position = "identity") +
  geom_vline(xintercept = optimal_threshold, color = "#D7191C", 
             linetype = "dashed", size = 1) +
  labs(
    title = "预测概率分布",
    subtitle = paste("最佳阈值 =", round(optimal_threshold, 3)),
    x = "预测为阳性(1)的概率",
    y = "频数",
    fill = "真实类别"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 11),
    legend.position = "bottom"
  ) +
  scale_fill_manual(values = c("#2C7BB6", "#D7191C"), 
                    labels = c("阴性(0)", "阳性(1)"))

print(p_prob)

# 将两个图形并排显示（可选）
grid.arrange(p_roc, p_prob, ncol = 2)
交叉验证
# 查看原始因子水平
cat("原始因子水平：", levels(Ttrain$`营养不良风险`), "\n")

# 方法1：自动转换为合法名称（如 "X0", "X1"）
# new_levels <- make.names(levels(Ttrain$`营养不良风险`), unique = TRUE)

# 方法2：手动指定有意义的名称（推荐）
# 假设原始水平顺序为 c("0", "1")，将"0"映射为"NoRisk"，"1"映射为"Risk"
original_levels <- levels(Ttrain$`营养不良风险`)
new_levels <- c("NoRisk", "Risk")   # 根据你的实际含义修改

# 重命名训练集的目标变量
Ttrain$`营养不良风险` <- factor(Ttrain$`营养不良风险`,
                                 levels = original_levels,
                                 labels = new_levels)

# 重命名测试集的目标变量（保持完全一致）
Ttest$`营养不良风险` <- factor(Ttest$`营养不良风险`,
                                levels = original_levels,
                                labels = new_levels)

cat("修改后的因子水平：", levels(Ttrain$`营养不良风险`), "\n")

# 设置随机种子
set.seed(123)

# 定义交叉验证控制参数（现在 classProbs = TRUE 可以正常工作）
ctrl <- trainControl(method = "cv",
                     number = 10,
                     savePredictions = "final",
                     classProbs = TRUE,
                     summaryFunction = twoClassSummary,
                     verboseIter = FALSE)

# 训练交叉验证模型
set.seed(123)
cv_model <- train(`营养不良风险` ~ .,
                  data = Ttrain,
                  method = "rf",
                  trControl = ctrl,
                  ntree = 200,
                  importance = TRUE,
                  metric = "ROC")
library(caret)          # confusionMatrix
library(pROC)           # roc, ci.auc, auc
library(ggplot2)        # ggplot绘图
library(ResourceSelection) # hoslem.test

# ================================
# 2. 准备训练集信息（假设已存在）
# ================================
# 训练集数据框 Ttrain 和随机森林模型 rfcla 必须已在全局环境中
train_levels <- levels(Ttrain$`营养不良风险`)   # 例如 c("0","1") 或 c("低风险","高风险")
train_features <- names(Ttrain)[names(Ttrain) != "营养不良风险"]

# ================================
# 3. 外部数据集（qwer）预处理
# ================================
# 因子水平对齐
if (!"营养不良风险" %in% names(qwer)) {
  warning("外部验证集 qwer 不包含 '营养不良风险' 标签列，只能进行预测，无法计算性能指标。")
  qwer$`营养不良风险` <- NA
} else {
  qwer$`营养不良风险` <- factor(qwer$`营养不良风险`, levels = train_levels)
}

# 检查特征列是否完整
missing_cols <- setdiff(train_features, names(qwer))
if (length(missing_cols) > 0) {
  stop("外部验证集缺少以下特征列：", paste(missing_cols, collapse = ", "))
}
qwer <- qwer[, c(train_features, "营养不良风险"), drop = FALSE]

# ================================
# 4. 模型预测
# ================================
qwer_pred <- predict(rfcla, qwer)           # 预测类别
qwer_prob <- predict(rfcla, qwer, type = "prob")  # 预测概率
cat("\n预测结果频数分布：\n")
print(table(qwer_pred))

# ================================
# 5. 性能评估（如果有真实标签）
# ================================
if (!all(is.na(qwer$`营养不良风险`))) {
  
  # 确保预测类别与真实标签水平一致
  qwer_pred <- factor(qwer_pred, levels = train_levels)
  
  # 混淆矩阵及性能指标
  cm_qwer <- confusionMatrix(qwer_pred, qwer$`营养不良风险`)
  cat("\n========== 外部验证集混淆矩阵 ==========\n")
  print(cm_qwer$table)
  cat("\n========== 外部验证集性能指标 ==========\n")
  cat(paste("准确率 (Accuracy):", round(cm_qwer$overall['Accuracy'], 4), "\n"))
  cat(paste("灵敏度 (Sensitivity):", round(cm_qwer$byClass['Sensitivity'], 4), "\n"))
  cat(paste("特异度 (Specificity):", round(cm_qwer$byClass['Specificity'], 4), "\n"))
  cat(paste("阳性预测值 (Precision):", round(cm_qwer$byClass['Pos Pred Value'], 4), "\n"))
  cat(paste("F1分数:", round(cm_qwer$byClass['F1'], 4), "\n"))
  
  # ROC曲线和AUC（假设阳性类别是 train_levels[2]）
  pos_label <- train_levels[2]   # 阳性类别（如 "1" 或 "高风险"）
  pos_index <- which(colnames(qwer_prob) == pos_label)
  if (length(pos_index) == 0) pos_index <- 2  # 备用方案：取第二列
  roc_qwer <- roc(qwer$`营养不良风险`, qwer_prob[, pos_index], 
                  levels = train_levels, direction = "<")
  auc_qwer <- auc(roc_qwer)
  auc_ci_qwer <- ci.auc(roc_qwer)
  cat(paste("\nAUC值:", round(auc_qwer, 4), "\n"))
  cat(paste("AUC 95% CI:", round(auc_ci_qwer[1], 4), "-", round(auc_ci_qwer[3], 4), "\n"))
  
  # 绘制ROC曲线
  roc_data_qwer <- data.frame(
    fpr = 1 - roc_qwer$specificities,
    tpr = roc_qwer$sensitivities
  )
  p_roc_qwer <- ggplot(roc_data_qwer, aes(x = fpr, y = tpr)) +
    geom_line(color = "#2C7BB6", size = 1.5) +
    geom_abline(linetype = "dashed", color = "gray50") +
    labs(title = "外部验证集 ROC 曲线",
         subtitle = paste0("AUC = ", round(auc_qwer, 4),
                           " [95% CI: ", round(auc_ci_qwer[1], 3), 
                           "-", round(auc_ci_qwer[3], 3), "]"),
         x = "假阳性率 (1 - 特异度)", y = "真阳性率 (灵敏度)") +
    coord_equal() + theme_minimal()
  print(p_roc_qwer)
  
  # ================================
  # 6. Hosmer-Lemeshow 校准度检验
  # ================================
  # 将真实标签转为0/1（1 = 阳性）
  y_obs <- ifelse(qwer$`营养不良风险` == pos_label, 1, 0)
  pred_prob_pos <- qwer_prob[, pos_index]
  
  # 动态调整分组数 g：避免因样本量太小导致负自由度
  n_obs_ext <- length(y_obs)
  g_max <- floor(n_obs_ext / 10)   # 每组至少10例
  g_use <- if (g_max >= 10) 10 else max(3, g_max)  # 至少3组，至多10组
  
  hl_test <- hoslem.test(y_obs, pred_prob_pos, g = g_use)
  cat("\n========== Hosmer-Lemeshow 校准度检验（外部验证集） ==========\n")
  cat(paste("分组数 (g):", g_use, "\n"))
  cat(paste("卡方值 (Chi-squared):", round(hl_test$statistic, 4), "\n"))
  cat(paste("自由度 (df):", hl_test$parameter, "\n"))
  cat(paste("P值:", format.pval(hl_test$p.value, digits = 4), "\n"))
  if (hl_test$p.value > 0.05) {
    cat("结论: P > 0.05，模型校准度良好（预测概率与实际观测无显著差异）\n")
  } else {
    cat("结论: P ≤ 0.05，模型校准度欠佳（预测概率与实际观测存在显著差异）\n")
  }
  
  # ================================
  # 7. 校准曲线 (Calibration plot)
  # ================================
  # 将预测概率分成等频组（与HL检验分组一致）
  cal_data <- data.frame(obs = y_obs, pred = pred_prob_pos)
  cal_data$group <- cut(cal_data$pred, 
                        breaks = unique(quantile(cal_data$pred, probs = seq(0, 1, length.out = g_use+1), na.rm = TRUE)),
                        include.lowest = TRUE, labels = FALSE)
  
  # 计算每组的平均预测概率和实际阳性率
  cal_summary <- aggregate(obs ~ group, data = cal_data, FUN = mean)
  cal_summary$pred_mean <- aggregate(pred ~ group, data = cal_data, FUN = mean)$pred
  
  p_cal <- ggplot(cal_summary, aes(x = pred_mean, y = obs)) +
    geom_point(size = 3, color = "darkred") +
    geom_abline(linetype = "dashed", color = "gray50") +
    labs(title = "外部验证集校准曲线",
         x = "预测概率 (均值)", y = "实际阳性率") +
    xlim(0,1) + ylim(0,1) +
    theme_minimal()
  print(p_cal)
  
} else {
  cat("\n外部验证集无真实标签，仅输出预测结果。\n")
  cat("预测为阳性(1)的概率统计：\n")
  print(summary(qwer_prob[, 2]))
}

# ================================
# 8. 与内部测试集性能对比（假设已有内部测试集结果）
# ================================
# 注意：内部测试集的 conf_matrix 和 auc_value 需提前计算并保存在环境中
if (exists("conf_matrix") && exists("auc_value")) {
  internal <- data.frame(
    数据集 = "内部测试集",
    准确率 = conf_matrix$overall["Accuracy"],
    灵敏度 = conf_matrix$byClass["Sensitivity"],
    特异度 = conf_matrix$byClass["Specificity"],
    AUC = auc_value
  )
  external <- data.frame(
    数据集 = "外部验证集",
    准确率 = cm_qwer$overall["Accuracy"],
    灵敏度 = cm_qwer$byClass["Sensitivity"],
    特异度 = cm_qwer$byClass["Specificity"],
    AUC = auc_qwer
  )
  compare_df <- rbind(internal, external)
  cat("\n========== 性能对比表 ==========\n")
  print(compare_df, digits = 4)
} else {
  cat("\n未找到内部测试集评估结果（conf_matrix, auc_value），跳过性能对比。\n")
}
library(pROC)
library(ggplot2)

# ---------- 1. 创建 ROC 对象（模拟数据） ----------
set.seed(123)
response_internal <- sample(c(0,1), 200, replace = TRUE)
predictor_internal <- rlogis(200, location = response_internal, scale = 0.6)

response_external <- sample(c(0,1), 150, replace = TRUE)
predictor_external <- rlogis(150, location = response_external, scale = 0.7)

roc_internal <- roc(response_internal, predictor_internal)
roc_external <- roc(response_external, predictor_external)

# ---------- 2. 计算 AUC 和 95% CI ----------
auc_internal <- auc(roc_internal)
ci_internal  <- ci.auc(roc_internal)

auc_external <- auc(roc_external)
ci_external  <- ci.auc(roc_external)

# ---------- 3. 放入列表并绘图 ----------
roc_list <- list(
  内部测试集 = roc_internal,
  外部验证集 = roc_external
)

p_combined <- ggroc(roc_list, legacy.axes = TRUE) +
  geom_abline(linetype = "dashed", color = "gray50") +
  scale_color_manual(values = c("内部测试集" = "#2C7BB6", "外部验证集" = "#D7191C")) +
  labs(
    title = "内部测试集与外部验证集 ROC 曲线对比",
    x = "假阳性率 (1 - 特异度)",
    y = "真阳性率 (灵敏度)",
    color = "数据集"
  ) +
  theme_minimal() +
  theme(legend.position = c(0.75, 0.25))

# 添加 AUC 文本
auc_text_internal <- paste0(
  "内部测试集 AUC = ", round(auc_internal, 3),
  " [95% CI: ", round(ci_internal[1], 3), "-", round(ci_internal[3], 3), "]"
)
auc_text_external <- paste0(
  "外部验证集 AUC = ", round(auc_external, 3),
  " [95% CI: ", round(ci_external[1], 3), "-", round(ci_external[3], 3), "]"
)

p_combined <- p_combined +
  annotate("text", x = 0.75, y = 0.25, label = auc_text_internal,
           color = "#2C7BB6", hjust = 0, size = 3.5) +
  annotate("text", x = 0.75, y = 0.15, label = auc_text_external,
           color = "#D7191C", hjust = 0, size = 3.5)

print(p_combined)

roc_list <- list(
  内部测试集 = roc_internal,
  外部验证集 = roc_external
)

p_combined <- ggroc(roc_list, legacy.axes = TRUE) +
  geom_abline(linetype = "dashed", color = "gray50") +
  scale_color_manual(values = c("内部测试集" = "#2C7BB6", "外部验证集" = "#D7191C")) +
  labs(
    title = "内部测试集与外部验证集 ROC 曲线对比",
    x = "假阳性率 (1 - 特异度)",
    y = "真阳性率 (灵敏度)",
    color = "数据集"
  ) +
  theme_minimal() +
  theme(legend.position = c(0.75, 0.25))

# 添加 AUC 文本
auc_text_internal <- paste0(
  "内部测试集 AUC = ", round(auc_internal, 3),
  " [95% CI: ", round(ci_internal[1], 3), "-", round(ci_internal[3], 3), "]"
)
auc_text_external <- paste0(
  "外部验证集 AUC = ", round(auc_external, 3),
  " [95% CI: ", round(ci_external[1], 3), "-", round(ci_external[3], 3), "]"
)

p_combined <- p_combined +
  annotate("text", x = 0.75, y = 0.25, label = auc_text_internal,
           color = "#2C7BB6", hjust = 0, size = 3.5) +
  annotate("text", x = 0.75, y = 0.15, label = auc_text_external,
           color = "#D7191C", hjust = 0, size = 3.5)

print(p_combined)

# 使用 ggroc 绘制两条曲线
p_combined <- ggroc(roc_list, legacy.axes = TRUE) +
  geom_abline(linetype = "dashed", color = "gray50") +
  scale_color_manual(values = c("内部测试集" = "#2C7BB6", "外部验证集" = "#D7191C")) +
  labs(
    title = "内部测试集与外部验证集 ROC 曲线对比",
    x = "假阳性率 (1 - 特异度)",
    y = "真阳性率 (灵敏度)",
    color = "数据集"
  ) +
  theme_minimal() +
  theme(legend.position = c(0.75, 0.25))

# 提取 AUC 和 95% CI（注意：您已有的 auc_ci 和 auc_ci_qwer 是长度为3的数值向量）
auc_text_internal <- paste0(
  "内部测试集 AUC = ", round(auc_value, 3),
  " [95% CI: ", round(auc_ci[1], 3), "–", round(auc_ci[3], 3), "]"
)
auc_text_external <- paste0(
  "外部验证集 AUC = ", round(auc_qwer, 3),
  " [95% CI: ", round(auc_ci_qwer[1], 3), "–", round(auc_ci_qwer[3], 3), "]"
)

# 在图上添加 AUC 文本（位置可微调）
p_combined <- p_combined +
  annotate("text", x = 0.75, y = 0.25, label = auc_text_internal,
           color = "#2C7BB6", hjust = 0, size = 3.5) +
  annotate("text", x = 0.75, y = 0.15, label = auc_text_external,
           color = "#D7191C", hjust = 0, size = 3.5)

# 显示图形
print(p_combined)
library(pROC)
library(ggplot2)

# 假设您已经定义好以下对象：
# roc_list, auc_value, auc_ci, auc_qwer, auc_ci_qwer

# ---------- 绘图 ----------
p_combined <- ggroc(roc_list, legacy.axes = TRUE, size = 1.2) +
  geom_abline(linetype = "dashed", color = "gray50", size = 0.8) +
  scale_color_manual(values = c("内部测试集" = "#2C7BB6", "外部验证集" = "#D7191C")) +
  labs(
    title = "内部测试集与外部验证集 ROC 曲线对比",
    x = "假阳性率 (1 - 特异度)",
    y = "真阳性率 (灵敏度)",
    color = "数据集"
  ) +
  theme_minimal() +
  # 图例放到左上角（x=0.15, y=0.85），避免与文字重叠
  theme(legend.position = c(0.15, 0.85),
        legend.background = element_rect(fill = "white", color = NA))  # 加白色背景更清晰

# 构建 AUC 文本（原样）
auc_text_internal <- paste0(
  "内部测试集 AUC = ", round(auc_value, 3),
  " [95% CI: ", round(auc_ci[1], 3), "–", round(auc_ci[3], 3), "]"
)
auc_text_external <- paste0(
  "外部验证集 AUC = ", round(auc_qwer, 3),
  " [95% CI: ", round(auc_ci_qwer[1], 3), "–", round(auc_ci_qwer[3], 3), "]"
)

# 添加 AUC 标注（放到右下角空白区域）
p_combined <- p_combined +
  annotate("text", x = 0.70, y = 0.20, label = auc_text_internal,
           color = "#2C7BB6", hjust = 0, size = 5.5) +
  annotate("text", x = 0.70, y = 0.08, label = auc_text_external,
           color = "#D7191C", hjust = 0, size = 5.5)

# 显示图形
print(p_combined)


internal_metrics <- data.frame(
  数据集 = "内部测试集",
  准确率 = round(conf_matrix$overall["Accuracy"], 4),
  准确率_下限 = round(conf_matrix$overall["AccuracyLower"], 4),
  准确率_上限 = round(conf_matrix$overall["AccuracyUpper"], 4),
  Kappa = round(conf_matrix$overall["Kappa"], 4),
  灵敏度 = round(conf_matrix$byClass["Sensitivity"], 4),
  特异度 = round(conf_matrix$byClass["Specificity"], 4),
  阳性预测值 = round(conf_matrix$byClass["Pos Pred Value"], 4),
  阴性预测值 = round(conf_matrix$byClass["Neg Pred Value"], 4),
  F1 = round(conf_matrix$byClass["F1"], 4),
  AUC = round(auc_value, 4),
  AUC_下限 = round(auc_ci[1], 4),
  AUC_上限 = round(auc_ci[3], 4),
  stringsAsFactors = FALSE
)

# ------------------------------
# 提取外部验证集指标
# ------------------------------
external_metrics <- data.frame(
  数据集 = "外部验证集",
  准确率 = round(cm_qwer$overall["Accuracy"], 4),
  准确率_下限 = round(cm_qwer$overall["AccuracyLower"], 4),
  准确率_上限 = round(cm_qwer$overall["AccuracyUpper"], 4),
  Kappa = round(cm_qwer$overall["Kappa"], 4),
  灵敏度 = round(cm_qwer$byClass["Sensitivity"], 4),
  特异度 = round(cm_qwer$byClass["Specificity"], 4),
  阳性预测值 = round(cm_qwer$byClass["Pos Pred Value"], 4),
  阴性预测值 = round(cm_qwer$byClass["Neg Pred Value"], 4),
  F1 = round(cm_qwer$byClass["F1"], 4),
  AUC = round(auc_qwer, 4),
  AUC_下限 = round(auc_ci_qwer[1], 4),
  AUC_上限 = round(auc_ci_qwer[3], 4),
  stringsAsFactors = FALSE
)

# ------------------------------
# 合并并打印对比表
# ------------------------------
compare_all <- rbind(internal_metrics, external_metrics)

cat("\n========== 内部测试集 vs 外部验证集 性能对比 ==========\n")
print(compare_all, row.names = FALSE)

最佳阈值分布图
# 确定阳性类别
pos_label <- levels(qwer$`营养不良风险`)[2]   # 例如 "Risk"
pos_index <- which(colnames(qwer_prob) == pos_label)
if (length(pos_index) == 0) pos_index <- 2     # 备选

# 构建绘图数据
prob_data_ext <- data.frame(
  True_Class = qwer$`营养不良风险`,
  Probability_1 = qwer_prob[, pos_index]
)

# 检查是否存在真实标签
if (!all(is.na(qwer$`营养不良风险`))) {
  
  # 阈值确定：优先使用内部测试集的最佳阈值
  if (exists("optimal_threshold")) {
    threshold_use <- optimal_threshold
  } else {
    # 若内部阈值未定义，则基于外部验证集计算（仅作演示）
    roc_ext <- roc(qwer$`营养不良风险`, qwer_prob[, pos_index],
                   levels = levels(qwer$`营养不良风险`), direction = "<")
    optimal_idx <- which.max(roc_ext$sensitivities + roc_ext$specificities - 1)
    threshold_use <- roc_ext$thresholds[optimal_idx]
    cat("警告：未找到内部最佳阈值，已使用外部验证集重新计算。\n")
  }
  
  # 提取外部验证准确率（若已计算）
  if (exists("cm_qwer") && !is.null(cm_qwer)) {
    acc_ext <- round(cm_qwer$overall["Accuracy"], 4)
    acc_text <- paste0("准确率 = ", acc_ext)
  } else {
    acc_text <- NULL
    cat("警告：未找到 cm_qwer 对象，无法显示准确率。\n")
  }
  
  # 绘制概率分布直方图（按真实类别填充）
  p_prob_ext <- ggplot(prob_data_ext, aes(x = Probability_1, fill = True_Class)) +
    geom_histogram(alpha = 0.7, bins = 30, position = "identity") +
    geom_vline(xintercept = threshold_use, color = "#D7191C", 
               linetype = "dashed", size = 1) +
    labs(
      title = "外部验证集预测概率分布",
      subtitle = paste0("最佳阈值 = ", round(threshold_use, 3),
                        if (!is.null(acc_text)) paste0("  |  ", acc_text) else ""),
      x = paste("预测为", pos_label, "的概率"),
      y = "频数",
      fill = "真实类别"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
      plot.subtitle = element_text(hjust = 0.5, size = 11),
      legend.position = "bottom"
    ) +
    scale_fill_manual(values = c("#2C7BB6", "#D7191C"), 
                      labels = levels(qwer$`营养不良风险`))
  
  print(p_prob_ext)
  
} else {
  # 若无真实标签，仅绘制预测概率的直方图（不分颜色）
  cat("外部验证集无真实标签，仅显示预测概率分布。\n")
  p_prob_ext <- ggplot(prob_data_ext, aes(x = Probability_1)) +
    geom_histogram(alpha = 0.7, bins = 30, fill = "steelblue") +
    labs(
      title = "外部验证集预测概率分布",
      x = paste("预测为", pos_label, "的概率"),
      y = "频数"
    ) +
    theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5))
  
  print(p_prob_ext)
}
